import pandas as pd
import numpy as np
import gurobipy as gp
from gurobipy import GRB
from shapely import wkt
import os
import geopandas as gpd
import folium
import contextily as cx
from IPython.core.interactiveshell import InteractiveShell
from shapely.geometry import Point
print(os.getenv('GRB_LICENSE_FILE'))
params = {
"WLSACCESSID": 'c819885e-2631-47d7-8108-8392e184ea0c',
"WLSSECRET": '534f1fd3-e9f5-4cd1-ac31-080d6502abf5',
"LICENSEID": 2398019,
}
env = gp.Env(params=params)
# Create the model within the Gurobi environment
model = gp.Model(env=env)
Building decarbonization is a critical component of the multiple transitions needed to combat climate change. The Inflation Reduction Act of 2022 provides unprecedented incentives for households to install home improvements for electrification and energy efficiency, like installing heat pumps instead of fossil-fuel powered heaters and weatherization/insulation (https://www.nytimes.com/interactive/2023/climate/tax-breaks-inflation-reduction-act.html)
These efforts have been aided by technological advances in the form of fine-grained building-level data, recently available at nationwide scale. One such dataset - developed by BlocPower - uses data on building type and equipment (single -family/multi-family, heating and cooling systems) sourced from tax assessment records. It then uses them as inputs for energy modeling algorithms/software developed by the Department of Energy. Crucially, this allows estimation of counterfactuals, i.e. the reduction in energy usage that would result from targeted retrofits. BlocPower’s own core operation takes building data for a city, uses their proprietary algorithm to estimate the ideal intervention, potential cost of a retrofit, and projected reduction in energy use (it also provides financing options).
Finance is the key variable in these ambitious transition plans. Many of the new policies provide incentives like tax credits and rebates, instead of direct transfers to finance retrofits. However, it is useful to get an estimate of the total cost required to achieve a specified reduction in energy use and emissions from buildings. (This is also useful where governments are spending money directly to decarbonize their own building portfolios). https://www.whitehouse.gov/briefing-room/statements-releases/2022/12/07/fact-sheet-biden-harris-administration-announces-first-ever-federal-building-performance-standard-catalyzes-american-innovation-to-lower-energy-costs-save-taxpayer-dollars-and-cut-emissions/
The set of possible interventions for each building can consist of:
Despite some differences in these particular interventions, the basic structure requires an upfront cost to be spent on a building, after which its energy usage will go down by some factor (i.e. there will be energy savings). For our current purpose, let weatherization be the only building intervention we consider. We abstract away the exact financing mechanisms included in the IRA and BIL (rebates and incentives), and assume that the local government has a fixed budget to allocate for weatherizing buildings. The objective is to maximize energy savings - or equivalently, reduce emissions - subject to the constraint that the total cost cannot be larger than the budget. For this, we need recommendations of which buildings to target.
Linear programming is a well-established mathematical tool used for policy purposes, starting from optimizing military logistics during WW2 to optimizing the energy source mixture (see appendix). This problem formulation fits well with a linear programming framework. Specifically, this is an example of a Mixed-Integer Programming (MIP) problem, as the solution must be an integer (we can't weatherize 3.6 buildings).
df = pd.read_csv(r"D:\Work\Georgetown\acad\mdi\summer_research\bp\geocoded_kentucky_zip_buildings.csv") # Reference any table in this project
test = df[['geocoded', 'building_id' , 'building_type' , 'heating_fuel_type', 'total_site_energy_GJ',
"address", "area_sq_ft", "energy_use_intensity"]]
test = test[test['heating_fuel_type'] != 'Unknown']
test = test.dropna(subset=['heating_fuel_type'])
#test = test.sample(frac=0.8).reset_index()
test['geocoded'] = test['geocoded'].str.replace('(','').str.replace(')','') # Remove parentheses
test[['lat', 'lon']] = test['geocoded'].str.split(',', expand=True)
test['lat'] = test['lat'].astype(float)
test['lon'] = test['lon'].astype(float)
# Convert latitude and longitude to a point
test['geometry'] = test.apply(lambda row: Point(row.lon, row.lat), axis=1)
Google Doc: https://docs.google.com/document/d/1JikyyMgS6zkjCQZMJQrtv8sU_ut1Uf0n8xkclPbHVg8/edit
Paper: https://www.sciencedirect.com/science/article/pii/S0306261922010510
$i \in T$: Index and set of potential buildings to weatherize.
$c_{i} \in \mathbb{R}^+$: The cost of weatherizing building $i$.
$e_{i} \in \mathbb{R}^+$: The energy use of building $i$.
$w_{i} \in [0, 1]$: This variable is equal to 1 if building $i$ is weatherized; and 0 otherwise.
The complication is that the cost $c_{i}$ and energy savings $e_{i}$ is different for each building, based on the chosen intervention and the building's own characteristics. Installing a heat pump or weatherizing in a gas-heating building of 6,000 sq. ft. would have a different cost from installing the same equipment in an oil-heating building of 3,000 sq. ft. Depending on local labor and material costs, even the exact same project on comparable buildings would have different costs in Wichita and Ithaca.
Industry experts have the following input:
Engie (Andrew Ludwig): 'Because of the wide variations, a heuristic approach is probably the best you can do.
BlocPower (Ankur Garg): 'The cost estimation is a process which requires local data, personnel hours and itself has an expense associated with it'
The original paper by Heleno et al (2022) seems to use a heuristic approach, by calculating cost and savings factors from the Weatherization Assistance Program (WAP) for different types of building archetypes. https://www.sciencedirect.com/science/article/pii/S0306261922010510
For this current example, we assume the following factors. So, weatherizing a gas-heated building reduces energy use by 4% and costs 2000 USD.
Multiply the above factors with the cost and energy columns to obtain 2 new columns per building, ei (energy savings after weatherization) and ci (cost of weatherization).
test = test[test['heating_fuel_type'].isin(['Oil','Gas'])]
test['energy_savings'] = test.apply(lambda row: row['total_site_energy_GJ'] - (0.96 * row['total_site_energy_GJ'])
if row['heating_fuel_type'] == 'Oil'
else (row['total_site_energy_GJ'] -
(0.98 * row['total_site_energy_GJ']) if row['heating_fuel_type'] == 'Gas'
else row['total_site_energy_GJ']), axis=1)
test['cost'] = test.apply(lambda row: 3000
if row['heating_fuel_type'] == 'Oil'
else (2000), axis=1)
test = test.reset_index()
# Create a new model
m = gp.Model("mip1", env=env)
# Add variables
w = m.addVars(test.index, vtype=GRB.BINARY, name="w")
# Set objective
m.setObjective(gp.quicksum(w[i]*test['energy_savings'].iloc[i] for i in test.index), GRB.MAXIMIZE)
# Add constraint: sum of w[i]*c_i <= budget
m.addConstr(gp.quicksum(w[i]*test['cost'].iloc[i] for i in test.index) <= 3000000)
# Optimize model
m.optimize()
# Check optimization status
selected_indices = [i for i in test.index if w[i].x > 0] # Store selected indices
# Add 'opt' column to the dataframe
test['opt'] = np.where(test.index.isin(selected_indices), 1, 0)
# Calculate and print final energy savings and cost of intervention
energy_savings = sum(w[i].x * test['energy_savings'].iloc[i] for i in selected_indices)
cost = sum(w[i].x * test['cost'].iloc[i] for i in selected_indices)
print("Energy savings achieved:", energy_savings)
print("Cost of intervention:", cost)
print('Buildings weatherized:', test[(test.opt==1)].shape[0] )
# gdf = df[['building_id', 'geocoded', 'area_sq_ft', 'address', 'energy_use_intensity']].copy()
# gdf = gdf.round(2)
# # merge to optimized dataframe
# gdf = gdf.merge(test, left_on='building_id', right_on='building_id', how='inner')
# Convert the DataFrame to a GeoDataFrame
gdf = gpd.GeoDataFrame(test, geometry='geometry', crs='epsg:4326')
# create retrofit column
gdf['Retrofit?'] = gdf['opt'].apply(lambda x: 'YES' if x == 1 else 'NO')
#round energy savings column
gdf['energy_savings'] = gdf['energy_savings'].round(3)
m1 = gdf.explore("opt", #tooltip=False, #categorical=True,
tooltip=["address", "area_sq_ft", "energy_use_intensity",
"building_type", "heating_fuel_type", "energy_savings", "cost", "Retrofit?"])
# Set `preferCanvas` to optimize performance of map
m1.options["preferCanvas"] = True
#set a tile layer - OSM
folium.TileLayer("OpenStreetMap").add_to(m1)
m1